!
!     1jigen kasyo hendo keisan
!    
!
!

implicit none

real::H(60),Z(60),B(60),QB(60),US(60),TS(60),DZ(60)
real::KP(60),IE(60)
real::PTIM(6)
data PTIM/0.,0.,1.,4.,12.,24./
!*****************************************************

real::G=9.8
real::S=1.65
real::D=5./1000.
real::RAMDA=0.4
real::Q=1000.
real::SN=0.02
integer::RL=3000
integer::DX=100
real::EPS=0.001
real::DT=10.
real::ETIM=86400.

!------------------------------------------------------

integer::NJ
integer::I,J
real::TIM=0.0,WF,H0
integer::ITM=1
NJ=RL/DX+1
PTIM(2)=DT
do I=3,6
   PTIM(I)=PTIM(I)*3600.
end do

!*******************************************************************

do J=NJ,1,-1
   KP(J)=(NJ-J)*DX
   B(J)=100.
   Z(J)=KP(J)/1000.
   if(KP(J)>1000..and.KP(J)<=1500.)then
      Z(J)=Z(J)+(KP(J)-1000.)/1000.
   endif
   if(KP(J)>1500..and.KP(J)<2000.)then
      Z(J)=Z(J)+0.5-(KP(J)-1500.)/1000.
   endif
end do
!write (*,*)'1'
!
!*****   WF   *****

call calWF(D,WF)

!********************************************************
!write (*,*)'2'
200 continue
if(TIM>ETIM)go to 900

!
!*****   hutouryu   *****

H0=3.02
call calH(H,Z,B,Q,SN,H0,G,DX,EPS,IE,NJ)

!
!*****   US    *****
do J=1 ,NJ
   US(J)=sqrt(G*H(J)*IE(J))
end do
!write(*,*)'3'
!
!*****   QB    *****

call calQB(NJ,S,G,D,US,QB,TS)

!
!*****   print out    *****
!if(TIM==PTIM(1).or.TIM==PTIM(2).or.TIM==PTIM(3).or.TIM==PTIM(4).or.TIM==PTIM(5).or.TIM==PTIM(6))then
do J=1,6
	if(TIM==PTIM(J))then
  	 call POUT(NJ,TIM,ITM,H,Z,B,DZ,US,TS,QB,KP,IE,DT,SN,Q,WF)
	endif
end do
!
!*****   DZ   *****

call calDZ(NJ,RAMDA,DT,DZ,Z,B,QB,DX)
TIM=TIM+DT
go to 200
!write (*,*)'4'
!
!*************************************************************

900 continue

stop
end

!
!****************************
!    cal of WF
!****************************

subroutine calWF(D,WF)
implicit none
real::D,WF,DD,A,DDD
DD=D*100.
A=.0036/1617.
DDD=DD**3.
WF=(sqrt(2./3.+A/DDD)-sqrt(A/DDD))*sqrt(1617.*DD)
WF=WF/100.
return 
end

!
!*****************************
!     hutouryu
!*****************************

subroutine calH(H,Z,B,Q,SN,H0,G,DX,EPS,IE,NJ)
implicit none
real::H(60),Z(60),B(60),Q,SN,H0,G,EPS,IE(60)
integer::DX,NJ
real::QQ,SNN,HH,BB,H3,FD,FU,FH,DFDH,DFDH1,DFDH2
integer::JJ,J
!write(*,*)H0,Q
QQ=Q*Q
SNN=SN*SN
H(NJ)=H0-Z(NJ)
do 100 J=NJ,1,-1
   JJ=J+1
if(J==NJ) go to 110
	HH=H(JJ)**2
	BB=B(JJ)**2
	H3=H(JJ)**(10./3.)
	FD=H(JJ)+Z(JJ)+QQ/(2.*G*BB*HH)+SNN*QQ*DX/(2*BB*H3)
	H(J)=H(JJ)
120 continue
	HH=H(J)**2	
	BB=B(J)**2
	H3=H(J)**(10./3.)
	FU=H(J)+Z(J)+QQ/(2.*G*BB*HH)-SNN*QQ*DX/(2*BB*H3)
	FH=FU-FD
	if(abs(FH)>EPS)then
   		DFDH1=QQ/(G*BB*H(J)**3)
  	 	DFDH2=(5./3.)*SNN*QQ*DX/(BB*H(J)**(13./3.))
  	 	DFDH=1.-DFDH1+DFDH2
   		H(J)=H(J)-FH/DFDH
   		go to 120
	endif
110 continue
IE(J)=SNN*QQ/(B(J)**2*H(J)**(10./3.))
100 continue
return 
end

!
!***********************************************
!     print out
!***********************************************

subroutine POUT(NJ,TIM,ITM,H,Z,B,DZ,US,TS,QB,KP,IE,DT,SN,Q,WF)
implicit none
integer::NJ
real::TIM,ITM
real::H(60),Z(60),B(60),DZ(60),US(60),TS(60),QB(60),KP(60),IE(60),DT,SN,Q,WF,PTIM(6)
integer::J
real::USPWF
character(len=80)::&
& fmt1 = '(10X A)',     &
& fmt2 = '(15X,A2,F10.1,A,5X,A2,F5.3,5X,A3,F10.1,A)',     &
& fmt3 = '(A5 F8.1 A1)',     &
& fmt4 = '(A5 F8.1 A2)',     &
& fmt5 = '(A,7X,A2,10X,A2,10X,A2,7X,A2)',     &
& fmt6 = '(I2,F6.0,F6.1,4F6.2,3E12.4,E12.4)'


write(*,fmt1)'*****  kasyo hendo keisan   *****'
write(*,fmt2)'Q=',Q,'(m3/s)','N=',SN,'DT=',DT,'s'
if(TIM<3600.)then
   write(*,fmt3)'TIME=',TIM,'s'
else
   write(*,fmt4)'TIME=',TIM/3600,'HR'
endif
write(*,fmt5)'NO KP B(m)   Z(m)   DZ(mm) H(m)  H+Z','IE','U*','T*','QB'
do 100 J=NJ,1,-1
   USPWF=US(J)/WF
   write(*,fmt6)NJ-J+1,KP(J),B(J),Z(J),DZ(J)*1000.,H(J),H(J)+Z(J),IE(J),US(J),TS(J),QB(J)

100 continue
ITM=ITM+1
900 continue

return 
end

!
!**************************************************
!     DZ
!**************************************************

subroutine calDZ(NJ,RAMDA,DT,DZ,Z,B,QB,DX)
implicit none
integer::NJ,DX
real::RAMDA,DT,DZ(60),Z(60),B(60),QB(60)
real::DQB,DZB
integer::JJ,J
do 100 J=1,NJ
   JJ=J-1
   if(J==1) go to 110
   DQB=QB(J)*B(J)-QB(JJ)*B(JJ)
   DZB=1./(1.-RAMDA)*DQB*DT/(DX*B(J))
110 continue
   DZ(J)=-DZB
   Z(J)=Z(J)+DZ(J)
100 continue
return
end

!
!******************************************
!     cal of QB
!******************************************
subroutine calQB(NJ,S,G,D,US,QB,TS)
implicit none
integer::NJ
real::S,G,D,US(60),QB(60),TS(60)
real::TSC=0.05
integer::J
do J=1,NJ
   TS(J)=US(J)**2/(S*G*D)
   QB(J)=8.*(TS(J)-TSC)**1.5*sqrt(S*G*D**3)
end do

return
end